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Firing patterns in the central nervous system often exhibit strong temporal irregularity and con¬ 
siderable heterogeneity in time averaged response properties. Previous studies suggested that these 
properties are outcome of the intrinsic chaotic dynamics of the neural circuits. Indeed, simplified 
rate-based neuronal networks with synaptic connections drawn from Gaussian distribution and sig¬ 
moidal non-linearity are known to exhibit chaotic dynamics when the synaptic gain (i.e., connection 
variance) is sufficiently large. In the limit of infinitely large network, there is a sharp transition from 
a fixed point to chaos, as the synaptic gain reaches a critical value. Near onset, chaotic fluctuations 
are slow, analogous to the ubiquitous, slow irregular fluctuations observed in the firing rates of many 
cortical circuits. However, the existence of a transition from fixed point to chaos in neuronal circuit 
models with more realistic architectures and firing dynamics has not been established. 

In this work we investigate rate based dynamics of neuronal circuits composed of several subpop¬ 
ulations with randomly diluted connections. Nonzero connections are either positive-for excitatory 
neurons, or negative for inhibitory ones, while single neuron output is strictly positive with out¬ 
put rates rising as a power law above threshold; in line with known constraints in many biological 
systems. Using Dynamic Mean Field Theory, we find the phase diagram depicting the regimes of 
stable fixed point, unstable dynamic and chaotic rate fluctuations. We focus on the latter and char¬ 
acterize the properties of systems near this transition. We show that dilute excitatory-inhibitory 
architectures exhibit the same onset to chaos as the single population with Gaussian connectivity. 
In these architectures, the large mean excitatory and inhibitory inputs dynamically balance each 
other, amplifying the effect of the residual fluctuations. Importantly, the existence of a transition to 
chaos and its critical properties depend on the shape of the single- neuron nonlinear input-output 
transfer function, near firing threshold. In particular, for nonlinear transfer functions with sharp 
rise near threshold, the transition to chaos disappears in the limit of a large network; instead, the 
system exhibits chaotic fluctuations even for small synaptic gain. 

Finally, we investigate transition to chaos in network models with spiking dynamics. We show that 
when synaptic time constants are slow relative to the mean inverse firing rates the network undergoes 
a transition from fast spiking fluctuations with constant rates to a state where the firing rates exhibit 
chaotic fluctuations, similar to the transition predicted by rate based dynamics. Systems with finite 
synaptic time constants and firing rates exhibit a smooth transition from a regime dominated by 
stationary firing rates, to a regime of slow rate fluctuations. This smooth crossover obeys scaling 
properties, similar to crossover phenomena in statistical mechanics. The theoretical results are 
supported by computer simulations of several neuronal architectures and dynamics. Consequences 
for cortical circuit dynamics are discussed. These results advance our understanding of the properties 
of intrinsic dynamics in realistic neuronal networks and their functional consequences. 


I. INTRODUCTION 


The firing patterns of circuits in the central nervous sys¬ 
tem often exhibit a high level of temporal irregularity. 
The effect can be seen by the inter-spike interval (ISI) 
distribution, which, except for a short refractory period, 
is similar to that of a Poisson process [1-3]. Intracellular 
recordings [4] indicate that this irregularity is due to fluc¬ 
tuations in the synaptic input to the neurons, suggesting 
a dynamic origin, and motivating the exploration of un¬ 
derlying neuronal circuit mechanisms. A possibly related 
issue is the ubiquitous diversity of the time averaged re¬ 
sponse properties of single neurons, e.g., their firing rates 
and tuning modulations, within a local population [1, 5]. 


Several theoretical studies explored the emergence of 
temporal irregularity and in particular chaotic dynamics 
in neuronal networks. These investigations focused on 
two types of models: rate-based models with Gaussian 
connections, and spiking dynamics of sparsely connected 
excitatory-inhibitory networks. The first class uses a fir¬ 
ing rate dynamics in which each unit is characterized by 
a smooth function that maps the synaptic input into an 
output firing rate. In its simplest version, the input- 
output transfer function is tanh(a:) where a zero value 
denotes some reference activity level and 1 and —1 de¬ 
note saturated firing and quiescent states, respectively. 
The architecture of the rate model was given by a ran¬ 
dom connectivity matrix where each connection is drawn 
from a Gaussian distribution, with zero mean and vari- 
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ance given by g'^/N, N being the size of the network. It 
was shown that the system exhibits a transition from a 
stable zero fixed point state for low value of g to chaotic 
state for large g. Furthermore, for large N this transi¬ 
tion is sharp and occurs at = 1 [61- In these models, 
the emergence of chaos is gradual, as both the ampli¬ 
tude of the fluctuations, their inverse time constant, and 
the Lyapunov exponent vanish as 5 —>■ 1“''. The chaotic 
state is asynchronous in that the correlations between 
fluctuations of different neurons are weak (and vanish as 
TV —>■ 00 ). 

The second class is motivated by biological reality. To 
capture the spiking dynamics, these models use either 
binary { 0 , 1 } neurons, or integrate-and-fire spiking neu¬ 
rons, and the connectivity is characterized by randomly 
sparse connections, where the mean number of connec¬ 
tions, K, is much smaller than N. To capture the biolog¬ 
ical constraints on the sign of the connections, the net¬ 
works consist of excitatory and inhibitory populations, 
where the nonzero output connections of excitatory (in¬ 
hibitory) neurons are positive (negative). The dynamics 
is dominated by the competition between strong excita¬ 
tory and inhibitory connections, leading to a dynamic 
cancellation of the excitatory and inhibitory inputs. The 
ensuing balanced state exhibits intrinsically generated 
Poisson-like stochasticity as well as asynchrony [7-10]. 
There is no parameter regime where the state can be 
characterized as a stable fixed point, and chaos does not 
emerge gradually as a function of synaptic strength. In¬ 
stead, temporal irregularity is always strong, and correla¬ 
tion times are short. The origin of the qualitative differ¬ 
ence in the behavior of the two types of models was never 
investigated comprehensively. In particular, it was un¬ 
clear whether the differences originate from the different 
dynamics (a smooth rate dynamics vs. spiking dynamics) 
or it is attributed to the different architectures: Gaus¬ 
sian distributed synapses in a single, fully connected, and 
statistically homogeneous population (with mixed excita¬ 
tion and inhibition) vs. two population (excitatory and 
inhibitory) architecture with sparse connectivity. 

The question of the existence of a bifurcation to chaos is 
not only interesting from dynamical systems perspective 
but may also have important functional consequences. 
Several studies have highlighted the computational util¬ 
ity of the nonlinear dynamics of random, recurrent net¬ 
works near the onset of chaos. For instance, a novel 
’reservoir computing’ model has been proposed, which 
utilizes the rich intrinsic network dynamics to learn to 
generate complex temporal trajectories [11-17]. Reser¬ 
voir computing is most effective above but near the tran¬ 
sition to chaos, due to the emergence of slow dynamical 
fluctuations, since many applications involve dynamics 
with time scales of seconds, much larger than the mi¬ 
croscopic time scales of a few milliseconds. Additionally, 
it has been shown that decoding signals from these net¬ 
works is particularly robust above and near the transi¬ 
tion to chaos [18, 19] . A recent study by Saxe et al [20] 
studies the dynamics of learning in deep networks. They 


define a critical point above which infinitely deep net¬ 
works exhibit chaotic percolating activity propagation, 
analogous to the chaotic state of recurrent networks. Fi¬ 
nally, recent advances in machine learning has generated 
resurgence of interest in recurrent networks (primarily of 
speech and language processing, e.g., [21-23] and refer¬ 
ences therein). Understanding the dynamics of generic 
recurrent networks will gain insight into these highly in¬ 
teresting computational capacities. 

In this work, we study the existence and the properties 
of the transition to chaotic dynamics in a broad range of 
models that span the two above mentioned model classes. 
In section II we introduce a general architecture for ran¬ 
dom recurrent networks with multiple sub-populations, 
obeying smooth rate based dynamics. We show the cor¬ 
respondence between randomly diluted networks in the 
balanced state and networks with Gaussian distributed 
connections with the same multiple population architec¬ 
ture. Section III introduces the mathematical framework 
of the dynamic mean field theory (DMFT) used in ana¬ 
lytical investigation of the properties of the network state 
and extend the theory of transition from fixed point to 
chaos, previously derived for a single gaussian population 
to the more general architecture. In section V we apply 
the theory to the simple case of a single inhibitory neu¬ 
ronal population and a threshold linear synaptic transfer 
function. The example of two population model is stud¬ 
ied in section VI. Although the DMFT is more complex 
than the single population network, we show that the two 
population network exhibits a transition from fixed point 
to chaos, that is similar to that of the single population 
case. The role of the single neuron nonlinear transfer 
function is elucidated in VII. First, we show that for suf¬ 
ficiently sharp function, chaos exists in large systems at 
all gain values. Secondly, the shape of this function de¬ 
termines the critical behavior of the relaxation times and 
largest Lyapunov exponent near the transition. 

Unlike the rate-based smooth dynamics, in models with 
spiking dynamics, a fixed point state does not exist as 
long as some of the neurons are firing. Nevertheless, it 
is interesting to explore the conditions under-which the 
spiking network exhibits a transition from a state with 
stationary inputs and firing rates, to a state where the 
underlying inputs and rates fluctuate in time. In sec¬ 
tion VIII we use a Poisson spiking model, to show that 
in the case where the synaptic integration time is much 
larger than the inverse of the neuronal firing rates, the 
synaptic current exhibits a sharp transition from a fixed 
point to chaotic dynamics, similar to that of a rate dy¬ 
namics. The behavior for large but finite synaptic inte¬ 
gration time is analyzed using scaling analysis, similar to 
that of a second order phase transition. The implications 
of the results for the understanding of the dynamics and 
computations in cortical circuits are discussed in section 
IX. 
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II. MODEL 

A. Randomly diluted network with multiple 
populations 

We consider a network of neurons composed of P sub¬ 
populations which are assumed for convenience to have 
equal size, N (Fig. lA). The recurrent connections, 
denote the synaptic efficacy between the presynaptic jth 
neuron of the Ith population to the post synaptic ith 
neuron of the fcth population, where k,l = 1 , ...,P, and 
The connectivity is randomly diluted, so 
that each connection is nonzero with probability p, 
where 

p = K/N (1) 

and zero otherwise. Thus, the mean number of inputs 
to each neuron is K. Some of the populations are as¬ 
sumed to be inhibitory. For an inhibitory population /, 
all are nonpositive. Conversely, all nonzero outgoing 
connections of an excitatory populations are positive. In 
addition, each neuron from the fc-th population receives 
a constant, uniform input, equal to nioWk, where the 
parameter niQ denotes the mean activity (firing rate) of 
neurons in the input population, and Wk are assumed 
positive. 

Throughout this work, we focus on networks with high 
degree of connectivity, i.e., N,K ^ 1. Although in most 
previous analytical work on the dynamics of dilute neu¬ 
ronal networks it was assumed that the network connec¬ 
tivity is sparse (i.e., 1 <C A' <C ) [ 8 , 10, 24, 25], here 
we will assume only that 1 K < N allowing for dense 
regime as well (in [9, 26] densely connected network were 
studied, but the focus was on the spatial patterns and 
correlations. The dynamical properties and the tempo¬ 
ral autocorrelations were not addressed). In order for the 
dynamics to be affected by the fluctuations in connectiv¬ 
ity, we assume all nonzero connections equal Jki/'/K, as 
discussed below. Similarly, the external connections scale 
as as Wk = uJk'/K- 



Figure 1. Network schematics. Arrows (circles) denote ex¬ 
citatory (inhibitory) connections. (A) The general network 
architecture studied here consists of multiple excitatory (E) 
and inhibitory (I) subpopulations driven by external inputs. 
Individual connections are randomly diluted or are Gaussian 
distributed. Subpopulations are distinguishable by the differ¬ 
ent statistics of connectivity. (B) The simplest model con¬ 
sists of a single population with random inhibitory recurrent 
connections with an external excitatory input. (C) Two in¬ 
terconnected randomly connected E and I populations. 


set all synaptic time constants to unity. The nonlinear 
transfer function ((>(h) denotes the firing rate of a neu¬ 
ron with an input synaptic current, h, and is analogous 
to the neuronal input-current to firing-rate transfer func¬ 
tion, known as the f-I curve. 


B. Rate dynamics 


C. Effective gaussian connectivity 


We first study a firing rate model for the dynamics of 
neuronal activity. The state of each neuron, say, with 
indices (i, k), at time t, is given by its total synaptic cur¬ 
rent that obeys a first order nonlinear dynamics, 

, p N 

= -hl{t) +^^JkA (^/ W) + (2) 

i=i j=i 

The first term is a decay term, the second is the synap¬ 
tic input from the network and the last is the synap¬ 
tic input from the external source. For convenience we 


Past work on random highly connected systems has 
shown that under fairly general conditions, in the limit of 
large K, the system’s behavior depends only on the first 
two moments of the connectivity matrix [27, 28]. Hence 
the connectivity matrix can be replaced by a ran¬ 
dom matrix of Gaussian distributed connections where 
the mean and variance are matched to that of the di¬ 
lute network. We thus consider a general dynamics of 
a multiple population network with fully connected con¬ 
nectivity matrix with Gaussian distributed connections. 
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where 


p N 


ki^ \'n 






1 

'mk{t) =—^(j){hl{t)) (4) 

' i=l 

are the mean population activities. The coefficients 
are quenched Gaussian variables with zero mean and vari¬ 
ance 


gki, and hence scales as y/K. On the other hand, the 
fluctuations scale as gki, hence are of order 1. Thus, it 
would seem that whenever the degree of connectivity, K, 
is high, the fluctuations induced by the random dilution 
will have negligible effect on the dynamics, reducing the 
system to that with a uniform connections and anoma¬ 
lously strong recurrent and external inputs. In fact, the 
system avoids this saturated state, by dynamically can¬ 
celing the mean excitatory contributions against the in¬ 
hibitory ones, resulting in the net mean inputs which are 
of order 1 , and of the same order as the fluctuations. 
Balanced states in networks with strong excitatory and 
inhibitory connections were previously studied (see [ 8 , 29] 
) in the context of binary or spiking networks. As we will 
show below, networks of units with smooth rate based dy¬ 
namics and transfer functions with finite gains also settle 
in this balanced state. 


{iJkif) = §- ( 5 ) 

We refer to g as the gain of the synaptic input. The 
contributions of the mean connections, is represented 
by the third term in the RHS of Eq (3) . 

In the dilute network, the mean connection between two 
populations equals pJki/'/K = '/KJ^i/N. The vari¬ 
ances of these connections are p{l — p)J^i/K = (1 — 
K/N)J‘^i/N. Hence the corresponding parameters of the 
equivalent Gaussian network are, 

9ki = iX-K/N)Jli, gki = '/KJm, hi = VKuJkmo. ( 6 ) 

Note that in the dilute networks, g^i and gki are related 
through 


Okl 9kl^ (7) 

In the theory below, we will analyze the dynamics of the 
Gaussian network, defined by Eq (3), in its generality; 
namely, we will consider the parameters describing the 
mean and the variance of the connections, gki and gki, as 
independent parameters. In addition, we will not restrict 
ourselves to the case that the mean connections gki and 
hi are large. The relation, (7), as well as the scaling of gki 
and hi with 0{\fK) will be adopted when applying the 
theory to the balanced dilute networks. The numerical 
simulations will use both Gaussian and random dilution 
and will show their equivalence. 


D. The balanced regime 

In the diluted network model, the net input from the 1-th 
population to each of the fc-th neurons, is proportional to 


III. FIXED POINTS AND THEIR STABILITY 

A. Mean field equations for the fixed points 

To fully characterize the state of the system we separate 
the population averaged quantities from the fluctuating 
ones, writing, 

h\{t) = Uk + 5hl(t) (8) 

where, 

+ hi (9) 

I 

The fluctuating inputs obey 

^Shlit) = -Shl{t) + Tfk{t). (10) 

where the spatiotemporal fluctuations in the currents, 
6hl.{t), are low pass temporal filters of the fluctuating 
synaptic inputs, 

=EE'^fc^''^(^i W) ■ (11) 

I j 

In the limit of a large number of inputs per neuron, these 
quantities obey Gaussian statistics with zero mean, zero 
spatial correlation but with temporal correlations which 
need to be evaluated self-consistently as explained below. 
First, we study the fixed point solution of the network 
dynamics, corresponding to a time independent state, 
where 5h\{t) = Shi = Sg], are Gaussian distributed with 




5 


zero mean and variance = {{Sh],)"^). This quantity is 
calculated self-consistently, as 


Afc = {{rfkf) -=^gliCu 


( 12 ) 


Note that due to the spatial summation, this quantity is 
essentially averaged over J. Interestingly, the response 
of the population averaged activity is coupled to the re¬ 
sponse of the population variances, defined as: 


where Ci = {(jP‘{h)) are the average autocorrelations of 
the neuronal activities. Finally using h\ = Uk + dril.^we 
obtain the self consistent equation for 


Ck = . (13) 

The constants Uk are determined self-consistently, via Eq 
(9), and 


Xki(.t) = 


mk = {(t> +M/c) y 


(14) 


Note that in Mean Field Theory all averages denote inte¬ 
gration over the Gaussian variables (z in the above equa¬ 
tions) which have zero means and unit variance. 


B. The balance equations 

In balanced architectures both gki and are of order 
Vk (see Eq (6)). In this case, the self consistent equa¬ 
tions for ruk assume a simple form. This is because for 
the system to settle into an unsaturated state, Uk must 
be of order 1; hence Eq (9) yield (to leading order in y/K) 


y^gkimi + hl=0, Vfc. 


(15) 


Since the mean rates are non-negative, Eq (15) can be 
obeyed only for a range of gki and values. Stability of 
the balanced state (see [8]) further restricts the parame¬ 
ter regimes. We refer to these restrictions on the param¬ 
eters as the balance conditions. Substituting the solution 
to the balance equations into Eq (14) yields equations for 
the residual, order 1 mean inputs, Uk- 


C. Stability of fixed points 

Stability of the population average activities The sta¬ 
bility equations for the population averaged degrees of 
freedom are determined by considering the response of 
the system to perturbations in the external fields, , 
which are uniform within the populations. It is conve¬ 
nient to define the uniform linear response by 


^ ^ ^ ^ dc^litydh^iO)- (16) 


dh°{0) 


Afc(t) = 


dhyO) 


{ShUt)6hl{t)). (17) 


(I - Ag -h iuj) 

-B 

' x(w)' 


A 

-Eg 

(I — D -I- iuj) 

X^{^\ 


E 


In the Appendix B we derive the following coupled equa¬ 
tions for the two sets of susceptibilities in the temporal 
Fourier domain. 


(18) 


Here, x y^are both PxP (P being the number of 
populations). The PxP dimensional matrices appearing 
in (18) are defined as. 


^ki — gki{4>i4>i) + {4>i) 

Bki = ^glM[My) + {Mi)] 

Dm= gh imr) + (Ml)) 

Em = 2gli{(j)i(j)'i) 

Thus, the fixed point is stable against population aver¬ 
age perturbations provided that all the eigenvalues of the 
matrix 


I-gA -B 
-gE I - D 


(19) 


have negative real part. 

Stability against local perturbations We now study the 
stability of the fixed point against small perturbations in 
the form of infinitesimal local fields . It is convenient 
to define the local susceptibility matrix 

X^i{t) = dhl{t)/dhfiO) (20) 

The average of the off-diagonal elements of this matrix 
is zero, and their variance is 0{1/N), hence we focus on 
the mean square susceptibility matrix, G, defined in the 
Fourier domain by 

Gfei(wi,,a;2) = N-^ Y^{Xm{^i)Xm{^2)) (21) 

■ij 

In Appendix C we show that the matrix G is 


G(a;i,,a; 2 ) = [(1 + + 1102 )! — M 


-1 


( 22 ) 
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where 

Mu = gli{i^'f). (23) 

Thus, the fixed points are stable provided the real parts 
of all eigenvalues of the P dim matrix M are all less than 
1. In general, M is not symmetric; however, the largest 
(in absolute value) of the eigenvalues is real. We will call 
M the stability matrix. 

In conclusion, we have derived two stability conditions: 
one related to population average perturbations, Eq (19), 
and a second related to local perturbations, associated 
with (22). Note that the mean interactions g do not 
appear in the latter condition. This is because the con¬ 
tribution to the ofl-diagonal elements of (21) from the 
uniform susceptibility is only of the order 1 /N . For the 
fixed point to be stable, both conditions must hold. How¬ 
ever the instability associated with each condition has a 
different implication. The population average instability 
signals either a runaway (as we will show in a concrete 
example in Section V) or a transition to another stable 
fixed point, a stable limit cycle, or some other coherent 
spatio-temporal states. On the other hand, as we will 
show below, the instability in (22) signals a transition an 
asynchronous chaotic state. 

Which of the instabilities occurs first when one varies 
one of the parameters depends on the specific architec¬ 
ture and parameter sets. Specific examples will be shown 
below. 


IV. CHAOTIC STATE: DYNAMIC MEAN 
FIELD THEORY 


Here, both y and 2 are Gaussian distributed random vari¬ 
ables with zero means and unit variances. The bound¬ 
ary conditions for the solution are: 9Afc(0)/9r = 0; 
9A^(oo)/9r^ = 0; and Afc(O) = A^ .The first condition 
stems from the general fact that A(t) is a symmetric, 
continuous function. The second one comes from the re¬ 
quirement that A(t) converges to a finite value as r —^ oo 
for a chaotic solution. Solution of (26) subject to these 
boundary conditions yields a unique solution for Akir) 
that converges to a fixed value, Afe(oo) at long time. Ex¬ 
cept for networks with h ^ —h symmetry, the quantity 
Afc(oo) is in general not zero, and represents the vari¬ 
ance of the Gaussian distribution of the time averaged 
synaptic currents. The fluctuations in these inputs de¬ 
termine the fluctuations in the time averaged firing rates 
of individual neurons. 

Details of the derivation of the DMFT are given in ap¬ 
pendix A. Alternative formalisms for deriving dynamic 
mean field equations for such architectures are Faugeras 
and coworkers using Mckean-Vaslov Fokker Plank for¬ 
malisms [30] , or for stochastic networks via path integrals 
(e.g., [31, 32]). 

Lyapunov exponent In the chaotic state, we expect the 
squared susceptibility to show an exponential divergence 
with time with a rate constant given by the largest Lya¬ 
punov exponent (LE). LE is defined as 

Xl= lim ^In^ (Gfci(T)), (27) 

r—>-oo ZT ^^ 
kl 

where 


The chaotic state is an asynchronous state with station¬ 
ary statistics, governed by the local-field autocorrelation 
functions 


Gm{t) 


lim 

dh°\t) J 


(28) 


A,{T)^{6hl{t)ShUt + T)), (24) 

where the angular brackets denote both average over 
neurons in the population and over time t. The self- 
consistent equations for these functions, derived from the 
Dynamic Mean Field Theory (DMFT) are (see Appendix 
A for detailed derivation): 

=J29kiGi{T), (25) 

where Cfe(r) = {4>(h],(t))(l)(h\(t + r))) are the firing 
rate autocorrelation functions, which depends on Akir) 
through 


Extending previous calculations [6] , we show in appendix 
G that the long time behavior of G(t) is determined by 
the ground state of a quantum mechanical problem with 
the Hamiltonian operator defined as 

^ = + (29) 

where 


Mkiir) =gli (^\Ja° - Ai(T)y + ^/A(t)z + 

(30) 

Finally, the LE is equal to 


Ckir) 



A, 


,(t)2/ -I- y^Ak{T)z + Uk^ 



Xl — —1 + Vl ~ eo) (31) 

where eo is the ground state energy of 'H. Note that, in 
general, "H is not Hermitian. However, complex values of 
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the Lyapunov exponents lead to unphysical oscillations of 
the squared susceptibility and so the ground state energy 
is expected to be real and negative. Direct rigorous proof 
is still missing. In the regime of a stable fixed point, M 
is time independent and eq > 0 , recovering the same 
stability criterion as above. A transition from a stable 
fixed point to chaotic dynamics occurs as eg vanishes. 


V. SINGLE INHIBITORY POPULATION WITH 
THRESHOLD- POWER LAW TRANSFER 
FUNCTION 

Solving the DMFT for general networks requires exten¬ 
sive numerics. Many of the salient features are captured 
by the simplest case of a single inhibitory population 
(Fig. IB) driven by a constant external excitatory in¬ 
put. The systems’s dynamics, 
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-h\t) = -h\t) + ^ {h^{t)) + -gm{t) + h°, (32) 

is characterized by the recurrent inhibitory mean gain 
parameter g < 0, the gain parameter of the synaptic 
fluctuations, g^ = and the excitatory exter¬ 

nal input, > 0. The system is further simplified by 
assuming a transfer function of a threshold- power law 
form [33, 34], 

^{x) = [x]’^, v>0, (33) 

where [a;]+ = max(a::, 0). For reasons explained below, we 
will see that the this monomial form is general enough 
when one studies the properties of the chaotic instability 
in networks with continuous threshold transfer functions. 


A. Fixed point and its stability 

In this model, the mean field equations for the fixed point 
are particularly simple. First, the mean input and its 
variance are given respectively, by 

u = gm + h^, (34) 

where m is the spatially averaged activity, and 

A = {6hj) = g^C, (35) 

where 5hi = hi — u and C is the mean square activity, 

m- 

The mean field equations for m and C are given by 

m = A'^/2^[z-ka;];), (36) 


and 

C = A"([z + x]^"). (37) 

Here 

X = u/'/A, (38) 

representing the mean input in units of the standard de¬ 
viation, and we have used the homogeneity property of 
the transfer function. 

Substituting (37) in Eq (35), the parameter x can be 
determined from 

l = f{[z + xY^), (39) 

where g = . The values of m and A are evalu¬ 

ated using Eqs (34), (38), and (36). 

As one increases g , the fixed point becomes unstable. 
Whether the first instability is the chaotic one, (Eq (22)) 
or the uniform one (Eq (18)), will generally depend on 
the parameter set (v, g and h^). 

Chaotic Stability A transition to a chaotic state occurs 
at the fixed point state when 

l = 5 V([z + a;]f-')). (40) 

For a threshold-linear transfer function {v = 1), g = g 
is independent of A and critical values gc and Xc, at 
which a bifurcation from a fixed point to chaotic dynam¬ 
ics occurs is determined by Eq (39) and (40) and are 
independent of h^. In this case C'{x) = H{—x) with 

H{x) = dze~^, and transition occurs when 

gc = v^iand Xc = 0. For non linear transfer functions, 
the critical values Xcand gc depend on the values of /i*^and 
9- 

Uniform Stability The stability condition for a network 
of single population against uniform perturbations is 
given by (See Appendix B where derive the uniform sta¬ 
bility condition in Eq (18) for a single population for a 
single population): 

+ + ( 41 ) 

Eor 1^=1, {(fxp") vanishes and the third term in the LHS 
of (41) is negative for g < 0, hence a chaotic instability 
(Eq 40) always occurs first (i.e., for lower values of g). 
For 1 / 1, the occurrence of uniform stability of the 

fixed point will depend on the parameters g and /i*^. 


B. Chaotic state 

In the threshold-linear model, it is useful to define the 
normalized autocorrelation 



g(T) = 1 - A(r)/A(0), (42) 

as well as 


g(oo) = 1 - A(oo)/A(0), (43) 

which measures the normalized variance of the dynamics, 
namely, the (spatially averaged) temporal variance of the 
local fields normalized by the squared amplitude. The 
autocorrelation function obeys (25) 


Stability of the chaotic solution The existence of a 
bounded chaotic phase, in which the mean activity of the 
network does not diverge and the trajectories of the local 
fields remain bounded requires stability of the uniform 
mode. Unfortunately a theory of uniform stability in a 
time dependent dynamical state is lacking. However, in 
the linear case, one can find simple arguments for the ex¬ 
istence of a chaotic solution. In this case, the normalized 
mean input a; < 0 in the chaotic phase is independent of 
h^, and one finds that only when 


l5l> 




(47) 


(44) 

Normalizing by Aq, and using the homogeneity of the 
transfer function as above, we can define a Newtonian 
equation of motion on the normalized autocorrelation 
(42) which reads 




dq ' 


(45) 


The potential V is given by 




I- -I 2 

oo oo 

9^ j V(1 - + a;) : 


— oo L-oo 


(46) 


where / Dx = f dx exp(x^ and $(/i) = 

fg dy<f{y) = [hYj^^ / {v-\-1) [6]. In the above equation we 

have used the equality Aq = A(0) and g = gA^ 

The initial conditions for (45) are g(0) = 0, dq{0)/dT = 0, 
and d‘^q{oo)/dT^ = 0. The second condition stems from 
the general fact that A(t) is a symmetric, continuous 
function. The last one comes from the requirement that 
A(r) converges to a finite value as r —>■ oo. These three 
conditions yield a unique solution and a unique value for 
the normalized mean input x. 

Due to the existence of a potential, x can be obtained 
without explicitly solving for qir) as follows. The above 
boundary conditions imply that the initial energy equals 
the potential energy, U(0), while the final energy equals 
the final potential energy, U(( 7 (oo)). Thus, conservation 
of total energy yields U(0) = U(g(oo)). Finally, for q(po) 
to be an equilibrium point, the force must vanish, hence 
dV{q{oo))ldT = 0. These two equations determine both 
q{oo) and x. Once x is known, (45) is integrated numer¬ 
ically to yield q{T). Finally, A(r) is evaluated by solving 
the mean field equation (36) with A = Aq. 


does a solution for the mean field equation exists; smaller 
values of |g|, entails dynamical instability. 

For non linear transfer functions, the above argument 
does not hold. For = 2 for example, a solution for 
the MF equation exist for any value of /I'^and g. How¬ 
ever, numerical simulations show that an instability in 
the chaotic phase exist, as can be seen in Fig. 2. 

We note that this instability results from the unbound¬ 
edness of (j). For (j) with a saturation level, the dynam¬ 
ics will always be bounded but, a crossover is expected 
from fluctuations spanning the linear dynamic range of 
the neurons when the net inhibition is large, to ’epilep¬ 
tic’ fluctuations in which neurons fluctuate between their 
saturated levels, for weak inhibition. 

Existence of a fixed point solution An interesting result 
that is implied by Eq (40) is that a stable FP exist only 
when 


> 1 / 2 . 


(48) 


In contrast, when i/ < 1/2, the RHS of the FP stability 
condition, (40), diverges indicating that no stable fixed 
point solution exists for finite g, and depending on g, the 
system is either in a stable chaotic state or diverges. This 
prediction is confirmed by the numerical simulations, Fig. 
2A, in which the normalized variance of the fields q{oo) is 
plotted as a function of For zz values of 0.5 or smaller 
the system is in a chaotic state even for small values of 
9- 

The instability of the fixed point for small v is due to the 
presence of positive local fields that are arbitrarily close 
to zero. Thus, for a system of finite size, where the posi¬ 
tive local fields are always of a nonzero minimum value, 
we expect that a fixed point will be stable at sufficiently 
small values of g. In the following sections, we focus on 
networks with threshold- linear (y = f) and quadratic 
(y = 2) transfer function, which exhibit a fixed point, a 
chaotic regime and a transition therefore. 

Finally we note that the same arguments hold for the 
multiple population case as well. 
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C. Phase diagram 


In Figure 2B we show the phase diagram for the 
threshold-linear {v = 1), depicting the regimes of sta¬ 
ble fixed point, chaos, and unstable dynamics in the pa¬ 
rameter space of g and g. For values of g < 1/2 the 
network settle into a fixed point. For larger gain values, 
if the inhibition is strong enough, i.e. the condition in 
(47) holds, then the dynamics is chaotic and bounded to 
a finite regime in the state phase. For lower values of 
the uniform inhibition, the mean activity diverges. Note 
that due to the semi-linearity of the transfer function, 
the phase diagram is independent of the magnitude of 
/iQ as long as it is positive. In a diluted network, the 
phase-plane is reduced to a single line 


9 = -9 


K 


1 - K/N ■ 


(49) 


For a sparse network where K ^ N, g = —'/K (dashed 
line in 2E). Thus, for large values of K, as in the balanced 
network case, the chaotic state at g > 1/2 is always stable. 
In 2C The phase diagram for the semi-quadratic ( 1 / = 2) 
is presented for = 1. For low values of uniform 
instability occurs at lower g than chaotic instability, and 
no stable chaotic phase exists. For larger values of 
a critical transition between a fixed point and chaotic 
dynamics exist. For larger values of gain g, the dynamics 
always diverges; this instability is shown by numerical 
simulation in the phase diagram. For a diluted network, 
where |g| = 0{'/K), there is always a chaotic transition, 
as shown in the inset of Fig 2C. 


D. Analytical evaluation of the Lyapunov exponent 

In the single population case, the evaluation of the Lya¬ 
punov exponent is also relatively simple. In this case, the 
single-component Hamiltonian is given by "H = + 

W{t) with a quantum potential W = —d^Vjdq^ where 
V is the classical potential (Fig 3A) , which can be easily 
evaluated once (/(r) is known (Fig. 3B). Indeed, for the 
semi-linear network we find that the ground state is neg¬ 
ative for g < 1/2 and positive otherwise, implying that 
Al is negative for g < y/2 and changes sign to a positive 
value for g > as expected in a chaotic state, see Fig. 
2C-F. For the non linear case within = 2 the Lyapunov 
exponents are calculated near and above transition in 
section VII . 


E. Numerical simulations 


Numerical integration of the full network equations in 
(32) verifies our theoretical predictions. For threshold- 
linear network, when g < y/2 the network settles into 





Figure 2. Phase diagrams for a single inhibitory population. 
(A) No stable fixed point for sharp transfer functions. Simu¬ 
lation results for the normalized variance q{oo) for networks 
with g — 0.2 and g = Q.l and different values of the power- 
law, V, characterizing the rise of the transfer function, Eq 
(33). Below V = 0.5 the variance is non-zero, implying there 
are temporal fluctuations even for low values of g. Simulation 
preformed on an inhibitory network [N = 5000) with ran¬ 
domly diluted connections and mean activity of m = 0.1. (B) 
Phase space diagram for a threshold-linear {u = 1) network 
of an inhibitory neurons with mean connectivity g < 0 and 
variance g . Dashed line shows an example for a transition in 
a diluted network (with K = 650 ) where because g = —y/Kg 
(see Eq (49)) the network lies on a line in the phase diagram. 

(C) Phase diagram for a threshold-quadratic {1/ = 2) network 
with an excitatory external field /i° = 1. Dashed line show 
uniform instability of the fixed point and solid line shows the 
transition to the chaotic phase (shaded area). Dotted (black) 
line shows the instability of the chaotic phase found by sim¬ 
ulations (Gaussian connectivity, N = 6000, h^ = 1). Inset 

(D) same as (C) for larger values of \g\. Dashed line shows 
existence line for a diluted network as in (B). 
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Figure 3. Calculation of the largest Lyapunov exponent for a 
threshold- linear network. (A) Numerical integration of equa¬ 
tion (45). The value of x is determined through the require¬ 
ment that the maximum of the potential at nonzero q equals 
its value at zero q. Compare the form of the potential with 
the exact value of x (solid line) with those calculated with x 
that deviates by ±1% from the correct value (dashed lines). 
(B) numerical solution for the normalized variance g(oo)as a 
function of g. (C)-(F) The normalized autocorrelation func¬ 
tion g(r) found by integrating the equation of motions (45) 
using the correct value of x for two values of g (top), and the 
corresponding quantum potential W = —d^Vjdq^ (bottom). 
The ground state energies eo (horizontal lines) were found nu¬ 
merically. In both cases they are negative implying a positive 
Lyapunov exponent. 


a fixed point, with the expected mean and variance of 
the local quenched fields (Fig. 4A). When g > \f2 chaos 
settles with temporal fluctuations that increase with g. 
The population averaged currents remain almost con¬ 
stant with small finite size fluctuations (Fig. 4B and 4C). 
The chaotic behavior is characterized first by the decay 
of the autocorrelation (Fig. 4D). which agrees with the 
theoretical (/(r), and second by a positive LE. The latter 
was calculated from simulations using Wolf’s algorithm 
[35]. The resultant values Al = 0.121, and 0.225 for 
g = 2.2 and 3.0, respectively, agree well with the values 
0.126 and 0.232 obtained by numerically calculating the 
ground state of the potential W{t) (Fig. 2D). Simula¬ 
tions of a semi-quadratic network also verifies our ana¬ 
lytical results. Simulations results near and above the 
chaotic transition are given below, in section VII. 

Randomly diluted networks The above numerical re¬ 
sults were for Gaussian distributed synaptic connections. 
However, as shown above, in the limit of large mean num¬ 
ber of connections per neuron, K, the DMFT is expected 
to hold also for randomly diluted networks, in which case. 


h t 




Figure 4. Numerical simulations of single population with lin¬ 
ear threshold transfer function. (A) Distribution of the local 
fields N in the fixed point regime. Thick curve shows nor¬ 
mal distribution given the theoretical mean and variance ob¬ 
tained from the mean field theory. Dashed (red) vertical line 
marks the firing threshold which is taken to be zero. (B),(C) 
activities of two networks with gains g = 2.1 and g = 2.8 
respectively. Bold lines show spatially averaged local fields; 
thin lines are local fields of a sample of four neurons. (D) the 
normalized AC function, A(r)/Ao = 1 — g(r) for two values 
of gain parameters as in (B) and (C). Solid black lines are 
the solutions of the equation of DMFT (45) superimposed on 
simulation results (dashed lines). Simulations preformed on 
a network (A = 6800) with Gaussian distributed connections 
and g = y/Kg, where K = 680 and external field hP = 1. 


g and g are related through (7) (Dashed lines in Fig 2). 
In the case of a single population with threshold linear 
transfer function, comparing the behaviors of different 
connectivity schemes (Gaussian, sparse and dense dilu¬ 
tion) is relatively simple since the normalized autocor¬ 
relation (/(r) depends only on g, not on g (see above). 
These expectations are borne out by our numerical sim¬ 
ulations, shown in Fig. 5. For simulations preformed on 
network with the same variance in their connectivity, the 
calculated normalized autocorrelation q^r) is identical in 
the Gaussian network and the randomly diluted networks 
with both p = K/N = 0.05 (sparse network) and p = 0.8 
(dense network). 


VI. TWO POPULATIONS WITH 
THRESHOLD-LINEAR TRANSFER FUNCTION 


Here we address briefly the application of the general 
theory to the case of a two-population network with one 
excitatory (denoted as E) and one inhibitory (/) popula¬ 
tions (Fig. 1C).. For simplicity we study the example of 
the semi-linear transfer function. Generalization to other 
values of i/ is straight forward, as presented in HI. 
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Figure 5. Robustness to changes in the synaptic distribution: 
(A) Autocorrelation function for sparse (p = 0.05), dense 
(p = 0.8) and gaussian connectivity Networks {N = 7000, 
h° = 1). The gains were chosen such that all networks have 
the same variance for the connectivity distribution {g = 2.8 
for the Gaussian distribution and g^/T^^ = 2.8 for the ran¬ 
domly diluted networks). (B) Normalized autocorrelation 
function 1 — g(r) = A(t)/Ao compared with the solution of 
the DMFT. 


A. Fixed point and its stability 


The fixed point equations for the variance and mean of 
the synaptic current are: 


Cfe('r) = V9k{T)y+ - qkiT)z + Xk 
for k G {E, /}, where 


J + 


(55) 


gfc(r) = 1 - Afc(r)/Afc(0) (56) 


are the normalized autocorrelation function for each pop¬ 
ulation, and here se set 


gki = gki 



Unlike the simple case of a single population, no classical 
potential function can be defined for the above two parti¬ 
cle motion in (54). Nevertheless, the dynamical equations 
can be integrated numerically, iteratively finding the val¬ 
ues for Xk that yield the desired asymptotic behavior of 
the normalized variance, 


Afc = ^ gliAi {[z + xi]^)^ , k = E,I, (50) 

l^EJ 

'mk = {[z + Xk] + )^, k = E,I, (51) 

alongside with the definition of Xk, 

Xk'/A^ = '^gkim + Wkirio. (52) 


qk{oo) = 1 - Afc(oo)/Afc(0). (57) 

Likewise, in general, the Hamiltonian governing the Lya¬ 
punov susceptibility G(t) is not Hermitian. However, 
as stated above, we expect the ground state to be real 
since the elements of G(r) are non negative by defini¬ 
tion, and complex ground state would mean oscillations 
around zero (See appendix C). 


C. Numerical results 


The eigenvalues of the stability matrix Mki = ghH^—xi) 
are given by 


A+ = 


Mee + Mil ± \J [Mee — Mji)^ + 4 :MieMei 


(53) 

Note that the two eigenvalues are real. For a fixed point 
to be stable A+ < 1. Thus, the transition to chaos oc¬ 
curs at parameters s.t. A+ = 1. When this eigenvalue 
becomes larger than 1, one must solve the DMF equa¬ 
tions for the chaotic state. 


Below we describe the chaotic state of this network, based 
on numerical investigations. In these simulations, we fo¬ 
cus on the balanced regime, in which all the mean con¬ 
tribution of each population to the local input is of order 
'/K. In this case, to leading order, excitation and inhi¬ 
bition cancel each other, and the mean activity of each 
population is given by 


ruE = -[J ^VFJbtoo (58) 

and 


B. Chaotic state 

The DMF equations can be written as, 

1 “ 9kECE{T) - ghCiir), (54) 


mi = —[J ^W]imo, (59) 

where W = [WeWi]’^ is the vector of feedforward con¬ 
nections, and J is the 2x2 matrix of the mean recurrent 
connections, as defined in the diluted model above. The 
balance conditions on the parameters are 
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det J < 0, 


JeiWi < JiiWe 


and 


JieWe < JeeWj. 

To demonstrate the properties of the two population sys¬ 
tem, we fix all the values of the connections except for a 
global gain g, and a, which controls the excitatory con¬ 
nections. For a = 0 the network is reduced to the single 
inhibitory network with the parameters used in section 
V. In Fig. 6 A we show a section of the phase space show¬ 
ing the critical gdo) line. The critical curve is calculated 
by solving the eigenvalues of the stability matrix (53) for 
each value of a. Fig. 6B, shows an example of the sta¬ 
bility for a = 0.55. Numerical simulations confirm the 
theoretical results as seen in Fig. 6C. For convenience of 
comparison, parameters were chosen so that they corre¬ 
spond to the same network parameters studied in [8, 9]. 
Unlike the binary network in which no stable fixed point 
exists for any g, in the threshold-linear network a transi¬ 
tion to chaos occurs at g = 1.21 (Fig 6B and C). Fig. 6C 
indicates that the normalized variance of the two popu¬ 
lations are very similar to each other (see discussion in 
Section 10 below). 

Figure 7 shows an example of the chaotic fluctuations in 
the two population network. Same parameters were used 
as in Fig. 6A and 6B with g = 1.6, within the chaotic 
phase.. As expected, inputs into neurons from both pop¬ 
ulations show large fluctuations, while the mean activity 
of each population is constant up to fluctuations of or¬ 
der Yj'/K (Fig. 7A). The autocorrelation of both pop¬ 
ulation decays monotonically with |r| to an equilibrium 
value (Fig. 7B). A signature of a dynamical balanced 
state is the substantial synchrony in the fluctuations of 
the excitatory and inhibitory mean activities. (Fig. 7C). 
Consistent with Fig 6C above, Fig. 7B shows that the au¬ 
tocorrelation functions of the two populations are nearly 
proportional to each other, implying that the normal¬ 
ized autocorrelation functions, ( 7 _e(t) and 9 /(t) are ap¬ 
proximately equal (as observed recently also by [36]; see 
Section VII below). 



Figure 6. Transition to chaos in a system of two popula¬ 
tions. parameters used: Jee = Jie = We = oig, Jii = 
—QjJei = —l.llf/ and Wi = 0.44y. Note that g is a global 
gain multiplying all the synapses; a denotes the excitatory 
efficacy. (A) The critical value of g as a function of the ex¬ 
citatory efficacy a. For 0 = 0 the network is identical to 
the single inhibitory network and phase transition occurs at 
g — y/2. (B) Largest eigenvalue of the stability matrix (53) 
for a network with a = 0.55 [marked by in panel (A)]. 
MF predicts a phase transition when A+ = 1 (asterisks). (C) 
Normalized variance, qk{oo), for inhibitory (blue triangles) 
and excitatory (red circles) populations as a function of the 
gain calculated from network simulations with a = 0.55. Cal¬ 
culated from simulation of a balanced diluted network with 
Ne = Ni = 3500, K = 700, connectivity parameters as in 
panel (B), and external activity mo = 1. Theory predicts a 
phase transition at g = 1.21 (asterisks), as seen in panel (B). 




VII. CRITICAL BEHAVIOR AT THE ONSET 
OF CHAOS 

In this section, we analyze the characteristics of the sys¬ 
tem at the onset of the chaotic state, and ask what fea¬ 
tures determine the critical properties of this transition. 
As will be seen below, these properties depend on the 
shape of the single neuron transfer function near the ori¬ 
gin. As such, the class of power law functions defined in 
33 can represent any continuous threshold function. 


Figure 7. (color) Fluctuations of an E-I balanced network with 
linear threshold transfer function. Same parameter set as in 
figure 6C. (A) Traces of the local fields of six neurons from ex¬ 
citatory (red) and inhibitory (blue) populations. Dashed lines 
show mean input into each population. (B) Time-lagged au¬ 
tocorrelation function of two population computed from the 
simulation. (C) Trace of the fluctuations in the spatially aver¬ 
aged activities 5mk{t) = mk{t) — {mkit)) ^ of both populations 
[same color codes as (A)[. Simulations of network similar the 
one used in Fig.6B, with g = 1.6. 
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Absence of stable chaotic phase Before we explore the 
critical properties near the transition to chaos, we note 
that not all transfer functions allow a stable chaotic so¬ 
lution. For example, Eq (25) may not have a solution 
satisfying all boundary conditions. In these instances dy¬ 
namic is either stationary or explosive (i.e., ’epileptic’). 
An example of such behavior is the exponential curve 
4>{x) = 6’”, as shown in Appendix D. On the other hand 
all 4 >{x) which saturate as a; —> oo are expected to exhibit 
a stable chaotic state for large values of g. 

In the following, we focus our study on transfer function 
of the type (33), with v > 1/2, which exhibits a phase 
transition from FP to stable chaos. 


A. Critical properties: A single population 





e Tj'fe 


When V > 1/2, the critical properties near the transi¬ 
tion to chaos depends on the value of z/. To see this, we 
first consider the case of a single inhibitory population. 
The critical value of g is given by the value at which 
g^iWiVAz + zr)]2) = 1, see Eq (39). As the chaotic 
state approaches the critical g from above, the ampli¬ 
tude of the time dependent fluctuations becomes small, 
hence, we can expand the MF equations in powers of ^(t) 
(which is small at all r). In addition, we will expand the 
equations in the small static parameter dx = x — Xc and 
the bifurcation parameter e = g'^ jg^-1. 

To leading order, the DMF equation (54) takes the form 
(see appendixE) 

= a{e)+ beq{T) + cqP{T), (60) 

where b and c are parameters of order unity, and the 
exponent p obtained from expansion of the firing rate 
autocorrelation, C(r) (see appendix E) is 


P = 





I 


(61) 


The constant term a vanishes at the transition and de¬ 
pends on e. In order for a nontrivial solution to exist, all 
the terms in (60) should be of the same order of e, hence 
the time scale of g(r) should scale as Tg// ~ l/-\/e- The 
amplitude of q{T) should scale as for i/ < 3/2 and as 
e for larger z/. Finally, a should scale as 0{eq) (as in¬ 
deed found by an explicit evaluation of a, see appendix), 
yielding the following scaling behavior, 


g(r) = ez.-!/, (62) 

where f{x) is of order 1. The Hamiltonian in (29) scales 
as 

W{t) = eF , (63) 


Figure 8. Critical behavior of a single inhibitory population 
with linear transfer function. Normalized variance (A), net¬ 
work relaxation time (B), and largest Lyapunov exponent 
(C) as a function of the distance from the critical point 
e = — 1. Circles show average over 20 simulations (Gaus¬ 

sian connectivity, N = 6000, h^ — 1) and dashed lines show 
theoretical predictions (See Appendix E). (D) Rescaled ID 
quantum potential F{t /Eq (63), for the Hamiltonian 
in (29). 


where F{x) is some other function of order 1. From (63) 
it implies that the largest Lyapunov exponent, Eq (31), 
scales as 


Xl = 0(e). (64) 


We have confirmed the above predictions numerically, for 
the cases of zz = 1 (Fig. 8 ) and v = 2 (Fig. 9). 

In this work, we assume a synaptic transfer function of 
the form (33). In [6] the variance of a network with a 
sigmoid synaptic transfer function, /)(/i) = tanh(h), and 

= 0, was shown to scale as A(r) = ef (“^^ and its LE 

as Al = O(e^). This behavior results from the h ^ —h 
symmetry of the system, which implies that Aq vanishes 
at the transition. This symmetry is not expected to hold 
in biological networks. Thus, the scaling shown in (62) 
and (64), in which the LE is larger than in the symmetric 
case, reflects the behavior in generic neuronal networks. 

In [37] the authors study the distribution of equilibria 
points above criticality, for a single population with sig¬ 
moidal transfer function. An interesting result is that 
the mean (with respect to realizations of Jy) number of 
equilibria behaves just like the Lyapunov exponent. Con¬ 
sequently, (64) may also elucidate the topological com¬ 
plexity of the flow above criticality for threshold power 
law transfer functions. 
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Figure 9. Critical behavior of a single inhibitory population 
with quadratic (v = 2) transfer function. Normalized vari¬ 
ance (A), network relaxation time (B) and LE (C) for small 
e = !gt — 1 gain above the critical point. Circles show aver¬ 

age over 20 simulations (Gaussian connectivity in the balance 
regime, N = 6000, K — 1200, /i° = 1) and dashed lines show 
theoretical predictions (See Appendix E). 


B. Multiple populations 


The above critical properties were derived for the case 
of a single population. However, we argue that these 
properties hold also for multiple populations. In the case 
of several populations, transition to chaos occurs when 
the largest eigenvalue of the stability matrix, M, (22), 
equals 1. Near the transition (on the chaotic side), this 
eigenvalue is slightly larger than 1 while the real part of 
the rest remain stable. To leading order, the unstable 
direction, , in the space of = 1 — Afc(T)/Afe(0) 

remains the same as the marginally stable eigenvector at 
the critical point. Hence, near the transition the domi¬ 
nant direction of temporal fluctuations are along the crit¬ 
ical eigenvector, and the dynamic equations collapse into 
one dimension, similar to the one population case. In 
general, we expect that all qk have nonzero components 
on the critical direction, hence gfe('r) Ri q{T)R^ and the 
critical properties are the same as the single population 
properties, above, see Appendix E. Similar arguments ap¬ 
ply to the scaling properties of the quantum Hamiltonian, 
hence the LE scales as (64). These results are supported 
by simulations of an excitatory - inhibitory network, de¬ 
fined in Section VI. The simulations displayed in Fig, 
10 demonstrate the universality of the critical proper¬ 
ties of the transition to chaos. In the specific case of a 
threshold-linear model, we find for a wide range of pa¬ 
rameters, R^ Ki (1,1,1, ...1), implying that near the tran¬ 
sition, the normalized autocorrelations of all populations 
are not only proportional but are expected to be approx¬ 
imately equal to each other, as detailed in appendix E. 
Note that the one-dimensional character of the chaotic 
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Figure 10. Critical behavior of a network with two populations 
and a linear transfer function. Diamond (red) and circles 
(blue) denote excitatory and inhibitory populations respec¬ 
tively. (A) The normalized variance of the two population, 
qE(oo) and qiipo), show the same critical behavior (a linear 
rise from zero) as the single population (Eq (62)). Due to the 
linearity, the normalized autocorrelations are nearly equal to 
each other near the critical point (see appendix E). (B) Net¬ 
work relaxation times of both populations are similar at the 
scaling limit. (C) LE of the entire network calculated over 
all populations, Eq (27), scale as a single population, Eq (64). 
Simulations conducted on network with diluted connectivity 
matrix {N = 3600, K = 600 for each population), and pa¬ 
rameter set as in Figs. 7 and 6C. 


fluctuations is exact only asymptotically close to the 
transition. Away from the transition, the non-linearity 
of C{q) couples the unstable mode to the other modes, 
inducing fluctuations in all modes. Interestingly, we find 
numerically that in many cases, even far from the tran¬ 
sition, the autocorrelations are still close to being equal, 
gfc(r) Ri qij) as is apparent in the example of Fig. 6C 
and 7B. Similar observations have been made recently in 
[36]. 


VIII. TRANSITION TO CHAOS IN SPIKING 
NETWORKS 


In this section, we inquire under what conditions the 
transition to chaos observed in the rate dynamics oc¬ 
curs also in spiking networks. In contrast to rate models, 
the membrane potentials of spiking networks are not at 
a fixed point as long as some of the neurons are active. 
Hence, in general, there is no transition from fixed point 
to chaos. An exception is the case of networks with slow 
synapses. If the synaptic time constant is long compared 
to the spiking dynamics, it is possible that for low synap¬ 
tic strength, the spiking dynamics is averaged out at the 
level of the slow synaptic currents, which will therefore 
stay approximately constant. Thus, there might be a 
critical g in which the state with almost constant synap- 
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tic currents undergoes an instability leading to a chaotic 
(or at least temporally irregular) state in which temporal 
fluctuations of the synaptic currents are large even at the 
scale of the synaptic time constant. 

To explore the possibility of transition to chaos in spiking 
networks with slow synaptic time constant, we consider 
a network of spiking neurons that fire with inhomoge¬ 
neous Poisson statistics, specified by the instantaneous 
rate of each neuron, r*(f) = cj) (^*(0)- We assume that 
the network is in an asynchronous state [38] and that 
the typical firing rate is large compared with the inverse 
synaptic time constant, 

Focusing on a single inhibitory populations, we write the 
dynamic equation of the local fields as, 

d ^ 

= -h\t) -bTogr(f) -b (65) 

i=i 

where is either Gaussian distributed connections or 
the corresponding randomly dilute network, as appears 
in Eq (32) above. The spike train {t) = ^ (i ~ 

of the presynaptic neuron j at times t\ represents the 
random Poisson process of neuron j. The non-linear rate 
function is given by the rectified linear transfer function 
(t>{x) = [x]_|_ where tq is a microscopic time constant, 

which is related to the inverse slope of the single neuron 
f — I curve. Note that in this model, the connections 
have units of time, hence we define 
so that g (and g) are dimensionless. For simplicity of 
notation, in the following we measure rates and times in 
units of To, i.e., we use units such that tq = 1. 

To understand the effect of spiking dynamics we can sep¬ 
arate the synaptic input into a rate contribution and a 
spiking noise contribution 

Vi{t) =V-{t)+g"P{t), ( 66 ) 

where 

N 

= ( 67 ) 

i=i 


N 



i=i 

(68) 

The auto-covariances of the two term on the RHS of 

are 

(66) 


(jllitWiit + T)) =g‘^C{T), 

(69) 

and 

{g7{t)ri7{t + T))=g\5{T), 

(70) 


where C'(r) = + t)) is the autocorrelation 

of the rate functions given by (26) (in units of Tq'^), 
and(5(T) is the Dirac delta function. The last equa¬ 
tion results from the average over the Poisson process, 

+ + = (j){hj{t))5{T) 

. Thus, the Poisson noise is equivalent to an additive 
white noise with amplitude g^r, where r is the popula¬ 
tion averaged rate r (in units of tJ”^). Thus, the spiking 
noise can be incorporated into the dynamic mean field 
theory, yielding, 

+ 5^^<5(r), (71) 

A. Perturbation analysis of the spiking noise 

As explained above, we are interested in the regime of 
large Ts (which is the synaptic time constant relative to 
To). This limit can be illuminated by writing the rescaled 
mean field dynamics, 

A(r)=g2C(f) + ^J(f), (72) 

where we have scaled time so that t = t/ts. The above 

equation demonstrates that the contribution of the Pois- 

2 

son noise (proportional to ^ ) is small relative to the 
rate contributions (which are proportional to g^r^ ) and 
the ratio between the two is of the order of (tTs)”^. 

For gain values above gc = and tts ^ 1, the noise 
contribution will be small compared to the unperturbed 
rate autocorrelation given by the solution of (45) (Figure 
llA). Below gc, the only time dependence comes from 
the spikes. To study this regime, we consider the spikes 
contribution as a small perturbation around the static 
autocorrelation, A(t) = Agt -b Asp{t), where Agt is the 
static solution for A in the rate model, and expand the 
equation above to linear order in Asp{t) yields 

(^1 - Asp{t) = g^C'A.p{t) + ^^(t), (73) 

where C = H{—x), see Eq (40). Solving (73) yields 

A.p(r) = ^e-^^, (74) 

where, following (40), 

X'^ = l-g^H{-x). (75) 

Here x is the unperturbed net input, Eq (38), and we 
have substituted back r = trs- [Note that g is in units 
of To[. Figure IIB displays the simulation results for the 
autocorrelation 6A{t) = A{t) — A(oo) for g below the 
critical value, and the theoretical estimates Asp(t). 
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Figure 11. Network of spiking neurons with inhomogeneous 
Poisson statistics. Simulation results for a network that fol¬ 
lows the dynamic equation in (65). (A) Above the critical 
point of the underlying rate dynamics, g^ = \/2, the rate 
term dominates the DMF equation, (72) and the autocor¬ 
relation calculated from simulation (rrs = 30 ,solid line) is 
fully explained by autocorrelations of the rates, found by (45) 
(dashed line). (B) For values of g well below the rate tran¬ 
sition, the rates do not have intrinsic fluctuations, and the 
fluctuations in the fields is solely due to the Poisson firing. 
MFT results of (74) (dashed line) agrees with simulations 
{rTs = 30, sold line). In intermediate values of g, rate and 
spike fluctuations interact nonlinearly. (C) Normalized vari¬ 
ance of the local fields q{oo) from simulations of networks at 
different gain levels in the transition region. For higher levels 
of rTs, the crossover between variance dominated by the Pois¬ 
son process and variance dominated by the fields dynamics 
becomes sharper. (D) The coefficient of variation {CV), Eq 
(76), for different values of gain. For low values of g, the CV 
approaches 1, indicating a Poisson process with (almost) con¬ 
stant rates. For higher values of the gain, the CV is larger, 
as expected from an inhomogeneous Poisson process. The 
crossover between the two regimes becomes sharper as rrs 
increases. Simulations preformed using an inhibitory pop¬ 
ulation with Gaussian distributed connectivity [N = 3500, 
K = 700) and g — hP — 1. 


B. Scaling analysis 


Approaching the transition point from below, the Poisson 
contribntion, (74), is amplified by the divergent response 
of the fields’ dynamics as A —^ 0. For a finite tTs how¬ 
ever, the autocorrelation remains finite at all g and the 
transition is smoothed. Indeed, Figure IIC shows that 
q{oo) increases smoothly as a function of g but its in¬ 
crease become sharper the larger rTgis is. Figure IID 
shows similar behavior for a measure of the fluctuations 
in the spike times, known as Coefficient of Variation, CV, 
defined as the ratio between the standard deviation of the 


ISI distribution, <jisi, and its mean, gisi [39],boundaty 

CV=^. (76) 

Pisi 

The CV increases from CV Ri 1, at small g, which the 
value for a homogeneous Poisson values, to substantially 
larger values above g = due to the fluctuations in 
the underlying rates. This increase is smooth but become 
sharper for large values of rr^. 

To study the effect of the small spiking noise on the tran¬ 
sition, one needs to perform a nonlinear analysis. Here 
we use a scaling analysis, similar to that of second or¬ 
der phase transition in equilibrium statistical mechan¬ 
ics. The scaling regime is characterize by two variables, 
rTs ^ 1 and e = g^/g^ — 1, |e| “C 1. We write the nor¬ 
malized variance near the transition as 

9(00) = 

where the scaling variable z is given by 


z = rTs\ef sign{e). (78) 

Far below the transition, z —?► —00 and according to (74) 
q~{oo) ~ requiring f{z -)■ -00) cx 

z°‘~^ and /3(1 — a) = ^. Similarly, above transition 2; —>■ 
00 , and from (62) q'^{oo) ^ entailing f{z —>■ 00 ) oc z“ 
and a/3 = 2. It follows that the critical exponents are 

4 5 

a = J,/S=5. (79) 


The behavior of the effective relaxation time of the net¬ 
work can also be treated in a similar manner. In the 
absence of the spiking noise, the effective time constant 
of the autocorrelation diverges as the transition is ap¬ 
proached from above, as Tefj ~ (Eq (62) and Fig. 

8B). In the presence of small spiking noise this time con¬ 
stant never diverges and we write. 


T-^ - 
^eff - 


(tTs) 


:F{z), 


(80) 


where Tg// is the effective correlation time in units of 
Ts- Here, the critical behavior both below and above 
transition, z —>■ ±cx), is similar (as can be seen from (74)) 
and ~ A ~ implying F[z —>■ ±(X)) (x \z\'^ and 

Pi = 5, or 



(81) 


Note that the above results predict that at the (rate dy¬ 
namic) transition (e = 0) the amplitude of the variance 
vanishes as (rrs)“^/® and the effective time constant di¬ 
verges {tTsY^^, respectively. Simulation results that sup¬ 
port these analytical predications are presented in Fig. 
12 . 
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Figure 12. Scaling behavior near criticality. Rescaled form 
for the normalized variance, {rTs)°‘q{oo), (A), and inverse 
effective relaxation time (rrs (B), using the scaling 

variable z = rsr\e\^sign{e). The data reveals the underly¬ 
ing scaling functions f{z) and F{z), respectively, consistent 
with (77) and (80) Dashed lines show the predicted asymp¬ 
totic behavior (see text). Simulations preformed on Gaus¬ 
sian network of various sizes {N = 2000, 3500 and 6000, 
h° = 1) and with various synaptic time scales (see legend). 
Effective relaxation was calculated from simulations by taking 
Teff = f T5A{T)dT/ f 5A{T)dr for r > 0. 


IX. DISCUSSION 

Universality of the transition to chaos across network ar¬ 
chitectures The results presented here show the univer¬ 
sality of the dynamical transition from a fixed point to 
chaos in large networks across different network archi¬ 
tectures. The theory, as well as simulations, show no 
dependence on the detailed distribution of the synaptic 
connectivity beyond the first two moments. In particu¬ 
lar, a Gaussian network behaves similarly to randomly 
dilute networks with either sparse or dense connectiv¬ 
ity. The two architectures differ in the mechanism by 
which fluctuations are not suppressed by the high con¬ 
nectivity. In the Gaussian networks, the mean of the 
connection distribution is assumed to be smaller by a 
factor of I/VN than their variance. In contrast, in the 
dilute networks, the mean inputs from each population 
is strong compared to the fluctuations and the fluctua¬ 
tions are nevertheless potent due to the dynamic balance 
of excitation and inhibition. This balancing amplifies not 
only the temporal fluctuations but also the time averaged 
ones. For instance, the neurons in the dilute networks ex¬ 
hibit a broad distribution of time averaged firing rates. 
In the threshold linear case, the distribution of the fir¬ 
ing rates will have a truncated Gaussian shape with a 
width which is related to Afc(oo). This balance mecha¬ 


nism is similar to the previously shown balanced state in 
spiking or binary networks. Here we show that this pro¬ 
cess is quite general and holds also in smooth dynamics 
with graded nonlinearity. It should be also noted that 
most of the previous work on balanced states assumed 
sparse connectivity, i.e., K <C N. In this limit, neu¬ 
rons share only few common sources, hence the network 
state is naturally asynchronous. Interestingly, we show 
here that even dense networks (where e.g. K/N is as 
large as 0.8, see Fig. 5) exhibit the same transition from a 
fixed point to asynchronous chaos, as in sparse or Gaus¬ 
sian networks (see also [9, 26]). It would be interesting 
to prove analytically the absence of stable synchronous 
states in these systems. 

The analytical and numerical investigations indicate that 
the transitions to chaos in a single population and a two 
populations (E-I) network is of the same nature and the 
DMFT predicts that this is true also for networks with a 
general multiple populations architecture. The extension 
of the DMFT to multiple populations is similar in spirit 
to previous work on stochastic multi-population networks 
with weak connections (of order l/K ) [40] and remains 
valid as long as the number of populations is small com¬ 
pared to K. The structure of the DMFT in the case of 
multi-population is different than in the single population 
case in that the dynamics of the autocorrelation function 
is not governed by a potential; thus, numerically solv¬ 
ing these equations is more challenging. Furthermore, 
although near the transition the autocorrelations behave 
qualitatively similar to the single population case, deep 
in the chaotic phase, depending on the parameter value, 
they may not be monotonically decreasing with time but 
might exhibit damped oscillations. Likewise, in multi¬ 
ple population networks, the ’quantum’ operator for the 
Lyapunov exponent is non Hermitian, the implications of 
which needs to be further explored. 

Stability against uniform perturbation Finally, we have 
focused primarily on asynchronous dynamical states and 
instabilities driven by the random components of the in¬ 
teractions. However, in the general architectures consid¬ 
ered here, the fixed point state may undergo first an in¬ 
stability associated with the population averaged pertur¬ 
bations. Interestingly, due to the inherent nonlinearity 
of our system, the uniform and local degrees of freedom 
are coupled. As a result, the response to uniform pertur¬ 
bations involves equations that couple the susceptibility 
of the population mean activities and that of the popula¬ 
tion activity variances, Eq.(18). This, together with (22) 
constitutes the stability conditions for the fixed points of 
our general architecture. These results can also be inter¬ 
preted as results regarding the eigenvalue spectrum of the 
random matrix of effective couplings which incorporate 
the local gains (ji'^h]) from linearizing around the fixed 
points, which are correlated with the random coupling 
(see also ]41]). 

In the examples studied here, the uniform instability 
leads to a runaway of the dynamics. In general, in the 
multi-population case, the fixed point instability to uni- 
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form perturbations may exhibit an instability leading to 
coherent oscillatory states, and these states may further 
become destabilized by the random connectivity yield¬ 
ing possibly chaotic states in the form of partially syn¬ 
chronous oscillations. The exploration of these states is 
beyond the scope of this paper. Finally, we emphasize 
that in the case of dilute networks, the primary focus of 
our work, the underlying strong inhibition puts the sys¬ 
tem aways from these instabilities as long as they are in 
states with excitation-inhibition balance. In particular, 
multiply stable states occur only under special conditions 
[29, 42]. 

Spatially modulated networks In the multi-population 
model, the functional structure of the entire network is 
given in terms of pairwise connectivity which depends 
on the subpopulation membership of the pair. In the 
brain, connectivity probability depends also on the dis¬ 
tance between neurons. This distance can be either with 
respect to physical location or functional, namely the pre¬ 
ferred stimulus of a neuron. Likewise, the external input 
to neurons is not in general homogeneous but may de¬ 
pend on the physical or functional coordinate of each 
neuron. Because these location dependencies are broad 
and smooth, the system can be well described by divid¬ 
ing it into columns, each represents many neurons having 
similar location. The total system, a hyper-column can 
thus be considered a special case of our general multi¬ 
population architecture. In the large N limit, the sums 
over columns will turn into smooth integrals over the spa¬ 
tial coordinates. In general, the architecture of a neural 
circuit if large enough will consist of both genuinely dis¬ 
crete subpopulations (e.g., Excitatory and Inhibitory) as 
well as continuously varying coordinates. An example 
is a model of orientation hyper-column in VI where the 
connectivity between neurons depend on the difference 
in their preferred orientation and the external input de¬ 
pends on the difference between the preferred orientation 
and the stimulus orientation [29]. In this case, the bal¬ 
ance equations (15) determine the tuning properties of 
the mean firing rates of the network. The DMFT de¬ 
scribes the statistical properties of the spatiotemporal 
fluctuations in the network dynamics. 

The shape of the nonlinear transfer function We have 
found that the existence of a transition from a fixed point 
to chaos as well as its critical properties when it exists, 
depend on the shape of the input-output transfer function 
near the origin (i.e., near the firing threshold). An impor¬ 
tant result is that for a transfer function rising as square 
root or sharper, i.e., for 4 ){x) oc x^'with v < 1/2, there is 
no stable fixed point and the system is chaotic also for 
small g. This raises the interesting question about the 
value of V in biologically relevant networks. In biological 
application of rate models, linear, quadratic or values of 
V larger than 1 have been often used [43-47]. The trans¬ 
fer function (p is often interpreted as reflecting the /-/ 
curve of a spiking neuron. The linear leaky integrate and 
fire (LLIF) model [48, 49] exhibits a sharp (1/| log(a:)|) 
rise from zero, corresponding to = 0. Our theory pre¬ 


dicts that random networks with such f-I curves exhibit 
chaotic dynamics also for low g. In conductance based 
(Hodgkin Huxley type) models the behavior of the f-I 
curve near threshold depends on the type and parameter 
values of the various nonlinear conductances contribut¬ 
ing to the spike generation. In Hodgkin Huxley (HH) 
Type I models, often used for modeling cortical neurons, 
the generic behavior of the f-I curve near threshold is a 
square root [50] making them border-line cases for a tran¬ 
sition at hnite g. Similar = 1/2 behavior is exhibited 
also by nonlinear (NLLIF) models [49]. The presence of 
adaptation current with long adaptation time constant 
results in a linearization of the f-I curve [51, 52], and 
thus corresponds to = 1. Thus, slow adaptation cur¬ 
rents in randomly connected networks are predicted to 
stabilize fixed point states at low g. It is interesting to 
note that the fixed point equations for the population 
averaged activities (14) are stable to linear perturbation 
even for i' < ^ due to the smearing of the singularity 
near threshold by the gaussian fluctuations. On the other 
hand, as we have shown here, this smearing is not strong 
enough to avoid instability of the fixed point of the local 
activities and fields. 

Finally, our prediction that for i/ < 1/2 no stable fixed 
point exists was derived in the limit of large K where 
the distribution of local inputs is Gaussian, hence un¬ 
bounded. However, in a finite system, where K is finite, 
there will be a finite gap between zero and the smallest 
input (in absolute value). Hence, in a finite system there 
should be a stable fixed point for sufficiently low values 
of g, even for low values of v. Studying this finite size 
effect is an interesting topic for future research. 
Although the properties of the transition to chaos are 
determined by the behavior of (j)(x) near threshold, the 
overall shape of (p may also affect the system’s behavior. 
For instance, in the threshold linear case, the effect of 
changing the magnitude of the external input is marginal, 
as it can be scaled out from the equations determining 
the chaotic behavior, due to the inherent linearity above 
threshold (see Eq (39)). In contrast, when <p{x) saturates 
to a finite value at large x, large external input pushes 
neurons to saturation and suppresses the onset of chaos 
[53]. Also, the unboundedness of the threshold-linear p 
leads to a divergent dynamics for sufficiently large g, see 
Fig. 2. Furthermore, when p{x) grows exponentially 
with X this instability sets in as soon as the fixed point is 
unstable. This divergent dynamics does not appear when 
p has a finite saturation value. 

Spiking dynamics The question whether networks with 
spiking dynamics exhibit a phase transition to chaos at fi¬ 
nite synaptic gain extends beyond the issue of the shape 
of the f-I curve. In contrast to recent claims [25], we 
have shown that a sharp transition from a fixed point to 
chaos in such networks is meaningful only when synap¬ 
tic time constant is large, where there is a clear sepa¬ 
ration between spiking dynamics and rate dynamics. In 
this limit, depending on the shape of the f — I curve, 
the underlying rates may be constant in low synaptic 
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gain and become chaotic above a critical gain, similar to 
the behavior of rate based dynamics. The general corre¬ 
spondence between smooth rate dynamic models and the 
dynamics of synaptic currents in neuronal systems with 
long synaptic time constants, has been studied previously 
[54] . However, the implications on the transition to chaos 
in random neuronal networks were not previously ad¬ 
dressed. In any realistic systems the time constants are 
finite, hence it is important to assess the smoothing of 
the transition due to finite synaptic times. Here we have 
characterized this smoothing by a scaling function with 
a single crossover exponent. This exponent determines 
the rate of change from stochastic spiking dynamics with 
static rates to smooth rate dynamics, as the synaptic 
time constant increases. This scaling analysis predicts 
relatively strong ’finite size’ effects of the spiking dynam¬ 
ics near the transition of the corresponding rate dynam¬ 
ics. First, the scaling regime is relatively large, given by 
\g — gc\ oc where r is the mean firing rate and 

Ts is the synaptic integration time. In addition, the ef¬ 
fective time constant of the synaptic fluctuations due to 
spiking, scales as Te// oc (tTsY^^. Thus, a sharp transi¬ 
tion requires rather large values of rrg. 

Concerning the biologically relevant regime, typical val¬ 
ues of mean rates range between an order of IHz to 
lOOHz. Fast synaptic currents (AMPA and GABAa) 
have decay time constants of the order of 1 — 10 msec so 
are not expected to exhibit the above transition. Slow 
synaptic currents (e.g., NMDA, GABAb, and peptide 
neurotransmitters) have time constants ranging from a 
few hundred milliseconds to minutes [39] and thus might 
be appropriate candidates for exhibiting transition from 
spike dynamics to rate fluctuations for systems with mod¬ 
erate rates (such as primate visual and motor cortex). 
Our analysis of networks of spiking neurons assumed 
inhomogeneous Poisson spiking model, which is a well 
known and extensively used statistical model of spiking 
fluctuations [48, 55]. However, it is interesting to ask 
whether our results hold for deterministic spiking mod¬ 
els with an appropriately smooth f — I curve. We believe 
our results regarding a transition to chaos for large K and 
large are valid also for conductance based spiking dy¬ 
namics, since they rely only on the separation of time 
scales between the firing dynamics and synaptic fluctu¬ 
ations. In particular, we expect that for weak synaptic 
gain the network will be in an asynchronous state char¬ 
acterize by synaptic inputs and population firing rates 
whose fluctuation amplitude is small on the scale of Ts. 
On the other hand, for strong synaptic gain, these fluc¬ 
tuations will be large even on the scale of Tg and the 
statistics of these fluctuation will follow the chaotic dy¬ 
namics of smooth rate dynamics. 

Beyond its simplicity, the advantage of the treatment the 
Poisson model is that it demonstrates the condition for 
rate chaotic fluctuations also for stochastically firing neu¬ 
rons. Our analysis of the Poisson model was restricted 
to threshold-linear rate function. It is straightforward to 
extend the DMFT equations to a general rate function, 


including one with a power-law behavior above threshold 
with ^ ^ 1 .It will be of interest to study in more detail 
the effect of fast spiking noise on the network dynamics, 
particularly for the cases where the (noiseless) f — I curve 
has < 1/2. 

One difference between the Poisson and the determinis¬ 
tic spiking dynamics is the degree of irregularity in the 
spike times, as measured e.g., by the standard deviation 
of the ISI distribution divided by its mean (known as the 
coefficient of variation-GV). In the Poisson model, the 
GV is close to 1 in the ’fixed point’ regime (as in homo¬ 
geneous Poisson process) and increases above it in the 
chaotic regime, due to the induced rate fluctuations. In 
contrast, in deterministic spiking models, the GV is ex¬ 
pected to be substantially lower than 1 in the ’fixed point’ 
regime. The detailed study of random networks with de¬ 
terministic spiking dynamics associated with sufficiently 
smooth f — I curve and slow synaptic time constants, is 
an interesting topic for future studies. 

In this work, we have addressed the effect of spiking dy¬ 
namics in smearing the transition from a state with static 
underlying synaptic currents to a state where the cur¬ 
rents themselves (hence the underlying rates) fluctuate 
in time. However, there can be other types of transi¬ 
tions from non-chaotic dynamics to chaotic dynamics in 
a spiking network even for fast synaptic time constants. 
In the case of a deterministic spiking networks, this tran¬ 
sition may mark the separation from a limit cycle at low 
g to a, chaotic attractor at high g. In this case, chaos is 
typically of fast time constant and cannot be described 
as instability in rates. In addition, it is likely that even 
in the non chaotic regime of deterministic spiking mod¬ 
els, irregular transients (with effective negative Lyapunov 
exponent, termed as stable chaos) will persist for a long 
time, and convergence time to the limit cycle will grow 
exponentially with the system size (as in [56-58]). Inter¬ 
estingly, the existence of chaos vs. stable chaos in spiking 
networks has been shown to depend on the details of the 
spike initiation dynamics [59]. A transition from irregu¬ 
lar but non-chaotic state to a truly chaotic state might 
be important for the information processing capacity of 
the system but will be hard to observe experimentally, 
as it is not reflected in the properties of the correlation 
functions. 

Our discussion of biological relevance of rate based dy¬ 
namics focused on identifying the units in the rate based 
models as single neurons and utilized temporal averag¬ 
ing as the mechanism for the adequacy of a rate based 
theory. An alternative scenario where rate description 
of a spiking network might be applicable is when the 
system consists of clusters of neurons, where spatial av¬ 
eraging can lead to rate dynamics (e.g. [52, 60, 61[). 

Under this interpretation, single units in our model rep¬ 
resent weakly synchronous neuronal subpopulations and 
the random connectivity corresponds to the large scale ef¬ 
fective connectivity between these populations. As such 
these models can serve as a useful framework for inves¬ 
tigating aspects of the large scale nonlinear dynamics of 
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the nervous system, as measured by EEG and fMRI sig¬ 
nals. 

This work explored the conditions under which random 
networks exhibit a transition from fixed point to chaos 
and the rate of development of chaotic fluctuations near 
the transition. These questions may bear important con¬ 
sequences for the computations that such network can 
produce. Recent models [13, 14, 62] have shown that 
random recurrent networks can be shaped through learn¬ 
ing to generate a broad repertoire of desired trajectories 
with biologically relevant time scales. These abilities pre¬ 
vail only near the onset of chaos. For substantially low 
gain, the activity in the network is strongly suppressed, 
whereas high above the transition, the intrinsic fluctua¬ 
tions are too fast and erratic to match the smooth and 
slow desired trajectories. It would be very interesting to 
study in detail how the results of the present work af¬ 
fect the computational powers of recurrent neuronal net¬ 
works. 

Note added: After the submission of this manuscript we 
became aware of a manuscript [63]that partially overlaps 
with the present work. 
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Appendix A: DMFT equations for the 
autocorrelations 

We extend previous derivations in ]6] to the current 
model of multiple subpopulation. Here, we show the 
derivation of the self-consistent dynamical equations for 
the population averaged autocorrelation function of the 
local fields, Afe(T), defined in (24). The dynamic equa¬ 
tions for the fluctuations in the local fields, 5h\.{t) , are 
given by (10), where in the large N limit, the noise terms, 
r7^(T), are Gaussian random variables with zero mean 
and {ql{t)r]l{t + T)) = J2i9kiCi{r), where Ckir) = 
{4>{h].{t))(j){h\,{t+T))) is calculated self-consistently by in¬ 
tegrating over a temporally colored Gaussian noise riKr). 
A convenient way of treating the temporal correlations 
between 6hk{t) and Shk{t + t) is to introduce three un¬ 
correlated Gaussian variables yi , y 2 and z with unit vari¬ 


ances such that 

5hk{t) = y/ayi -I- 

5hk{t + t) = y/ay2 + \fPz. 

with 

a = Afc(O) - Afe(T), 

/? = Afc(r). 

With these variables, the population averaged autocor¬ 
relation function can be written as the integral over in¬ 
dependent Gaussian variables 

Ckir) = 

J Dyi J Dy 2 j Dzcj) (^y/ayi + y/^z + Uk'^ x 

4> [yay2 + y/Pz J- . (Al) 

Finally, to yield the self-consistent equations for Aj, it is 
convenient to use the Fourier transform of the dynamic 
equation (10) one gets 

(1 -k iu;)Shl{uj) = yliuj). (A2) 

from which we obtain, 

(1 -f (\Shliu})\^'j = =J29kiCiiu:). 

(A3) 

Finally, preforming a Fourier transform back to the time 
domain, and substituting the result of (AI), the P self 
consistent equations for the autocorrelations of the local 
fields read 



'^9ki (\/Aj(0) - Ai(r)y -k y/ A,(t) 2; -k ui'j)■ 


Appendix B: Stability Equations for the Fixed 
Points 

In this Appendix we derive the stability conditions for 
the fixed point states of the network. 


1. Population average linear response 

For simplicity we present the derivation with one popula¬ 
tion architecture. The extension to multiple population 
is straightforward. 

We evaluate the time dependent response function 
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with t' —> oo. Finally, we note that in this limit, 


X^{t) = d^{h,{t))/dh°{0) (Bl) dA{t2t) d{Sh{t2)Sh{t)) 

dhO{0) dhO{0) 

and denote by x(t) the spatial average of Xi- 

at the fixed point solution. From the equations of motion. We denote 
we obtain 


(B7) 


(l+ 5 i)XiW = (l}'ihi)Y^JijXj{t)+(l)'{hi)gx{t)+(l)'{hi)S{t) 

J where, 

(B2) 

where (j)'{hi) are the derivatives of the activation function 
evaluated at the fixed points. Let us consider the general 
term 


= 2x^(oO)0 (B8) 


X^(t',0 


dA{t', t) 
dW{Q) 


(B9) 


</>'(/!.) E '^bX.W = E W) (B3) 

3 ^ 3 

Because the fixed point values of hi are independent of a 
perturbation at time 0 . Thus, using Eqs (10) and (11), 
we have 

(/>'(hz)E'^yWW = -Q^il + dt)(i)'{hi)6hi{t) (B4) 

where dt stands for the time derivative operator. After 
averaging and using mean field theory. 


{()>' {hi{t2))5hi{ti)) = 

A(t2, t2) - A{t2, ti)x + yjA{t2,tl)z + u{t2))x 
[i/A{ti,ti) - A{t2,ti)y + A{t2,ti)z)) 

= {4>'{VA{t 2 M) - A{t 2 ,ti)x+ 
i/A{t 2 ,ti)z + u{t 2 ))\/ A{t 2 ,ti)z) 

= A{t2ti){(j)”{\J A{t2,t2)z + u{t2))) (B5) 

where A{t,t') = {Sh{t)6h{t')). Thus, 

(0'(h,)E j',,x,(t)) = (i + ao^|^(<(>") (B6) 

1=1 ' I 


Note that the factor of 2 was introduced, because in the 
fixed point. 


X^{t,t) = 2x^{oo,t) = x^{t) (BIO) 

which is the susceptibility of the mean equal time vari¬ 
ance. 

Substituting the above results in (B2) and averaging we 
obtain, 

{l-{cj3')g+dt)x{t) = \{r){l+dt)x^{t)+{W) (Bll) 


Calculating of 

From the DMFT we obtain 


{l + dt,){l + dt,)A{t2ti) = 9^C!{ti,t2) (B12) 

where. 


<^(* 1 ,^ 2 ) = {4>{\/A{t 2 ,t 2 ) - A{t 2 ,tl)x+ 

y/Ait 2 ,ti)z + u{t 2 ))(l){\/A{ti,ti) - A{t 2 ,ti)y+ 

A{t 2 ,ti)z + u{ti))) (B13) 


Thus, 


(1 -f (9tJ(l -f dt^)x^{t2ti) = 9‘^[dC!{ti,t2)/dA{t2,ti)x^{t2,ti) + dC{ti,t2)/dA{t2,t2)x^{t2,t2)+ 

+ dC{ti,t2)/dA{ti,ti)x^{ti,ti) + {(I)'{t2)<i)iti))[gx{t2) + S{t2)] + {(t>{t2)(l)'{ti))[gx{ti) + S{ti)]] (B14) 


where, by derivation of (B13), 

dC{ti,t2)/dA{t2,ti) = {(j3'{t2)(l)\ti)) (B15) dC{h,t2)/dA{ti,ti) = ^{(t>{ti)(l)"{ti)) (B16) 
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and 


du 

Mo 


gx{t) + S{t) 


(B17) 


We will assume t 2 > ti > to are all large but their differ¬ 
ence is of order 1. 

Let us take t 2 — ti ^ oo . Then, substituting (BIO) one 
obtains: 


Multiple population network A straightforward general¬ 
ization to the multiple population case yields for 


and 


Xki{t) = 


dmkjt) 

dhi{0) 


(B27) 


Xkiit) 


d 

dhi{0) 


Afc(i). 


(B28) 


{^ + dt)xti{t) = 


Eqs (Bll) and (B18) constitutes two coupled equations 
for the response functions of the mean and variance of the 
population activity. It is perhaps convenient to eliminate 
the time derivative in the RHS of (Bll) by substituting 
in it Eq (B18), resulting in 


(1 - {ct>')g + dt)x{t) = + (#")]X^(0+ 

{^")g^{^mx{t) + S{t)] + {^')S{t) (B19) 


glm [{^m) + {4'm4’m)] Xmli'^) 

m 

+ 2 'y ^ gkm{4’m4>m) ^ ^ gmm'Xm'l {t) -f 5ml5{t) 

m \_m' 

(B29) 

or, in matrix notation, 

(I-D-fat)x^=E(gx + I(5(t)). (B30) 

Likewise, the dynamical equation on y(t) is 


We can write these equations by 

{a + dt)xit) = bx{t)+c6{t) 

{d + dt)x^{t) = exit) + fSit) 
or, in Fourier and matrix representation. 


(1 — ga -I- iiS) 

-b 

X 


a 

-ge 

{1 — d + iuj) 



e 


where 

a = (<^'V(#') + (<^'), 

and 

e = 2 g^{(j)(j)')g. 


(B20) 

(B21) 


(B22) 


(B23) 


(1 + dt)xkiit) = 

'^gkmi<P'm)Xmlit) + (())fe)4i <5(0 

m 

+ l{€)i^ + dt)x^,it) (B31) 
or 

il-Ag + dt)xit)=Bx^ + ASit). (B32) 

Where in the above we have defined the P x P matrices 
Afc/ = Skiifj})) + {4'k)gki{4>i4’'i)y (B33) 

Bki = ■^glii^'k) [('('f) + {4>i4>'i)] , (B34) 

Dki = gh M) + iMi)] , (B35) 

and 

Eki = 2glMp^i)- (B36) 

Eq (B30) and (B32) can be written using a 2P x 2P 
matrix in Fourier space as 


(B24) 


(B25) 


(B26) 


(I - Ag -1- iuj) 

-B 

x(w) 


A 

-Eg 

(I — D -I- iio) 

[x^(^)J 


E 


(B37) 


Thus, the fixed point is stable against population aver¬ 
age perturbations provided that all the eigenvalues of the 
matrix 


I-Ag -B 

-Eg 1-D 


(B38) 


have negative real part. 
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2. Stability of fixed point to local perturbations 

The susceptibility matrix of the network to a random 
local perturbation is defied as 


where h^\t') are infinitesimal external perturbations 
around the fixed point state. We are interested in the 
mean square of this quantity. Using the dynamical equa¬ 
tion (2), we find the time evolution of a the susceptibility 


Xkiit -t')= dhl{t)ldh^i\t') (B39) 

I 


P N 




E E (M-) *;.'/(< - f) + y Es«' E - t') + S{t - t')SijSkh ^ > 0 (B40) 




where Sij and S{t) are the Kronecker and Dirac delta 
functions respectively. Defining, 


Gfci(wi,,a;2) = iV ^ (B41) 

u' 

and averaging (B40) yields 


G(wi,,W 2) = [(l + iwi)(l + iw2)/-M]-i (B42) 

where M is the stability matrix, M^i = {4)\h].)(j)'{h\)) 
(Eq (23)). 

Note that Gki are of order 1; the contribution of the 
uniform Xki to them is of order 1/iV hence negligible. 
Furthermore, the uniform components of the interactions 
(in Eq. (B40)) do not contribute to (B41) to leading 
order as they are smaller by a factor of 1/iV relative to 
the random contributions. 


lution of a the susceptibility 



p N 

EE'7»U'(4'(<))y.'/(i.i')+ 

X 9ki' (^i' (*)) ^')+ 

i'=i j' 

S{t - t')SijSki, t > 0 (C2) 

where Sij and S{t) are the Kronecker and Dirac delta 
functions respectively, denoting a local spatiotemporal 
perturbation. 

Population averaged susceptibility Let us definey^ 
as the within-population spatial average susceptibility 

xhit = Xkiit^t') = (xl]it^t'))j 

i,j=l 

(C3) 

This uniform matrix is formally given by 


Appendix C: Stability Equations and the largest 
Lyapunov Exponent (LE) 


The susceptibility matrix of the network to a random 
local perturbation is defied as 


Xki^'^P) = dhl(t)/dh°\t') (Cl) 

where h^\t') are infinitesimal external perturbations. 
Using the dynamical equation (2), we find the time evo- 


x°it-t') = { 


1 -1 


/ _ 7_/O 


)a{t)S{t — t'), t>t' 
(C4) 

where / is spatial identitiy matris, a{t) is the tempo¬ 
ral operator a = {1 + d/dt)~^, (^/') 

•^kl^ ~ ^9ki<P' ■ [The RHS is formally an NPxNP 

matrix; however, it is meant to be reduced to a PxP 
matrix by averaging over the i,j indices]. 


Expanding the RHS in powers of J shows that all con¬ 
tributions vanish upon averaging, hence 
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)a{t)5{t — t'), t>t' 


(C5) 



xlM = ^9kk' {4>' {h'-i))xl'i{t) + ^5ki5{t), t > 0 

k' 

(C8) 


Furthermore, performing the averaging in the RHS of C5 
and fourier transforming the result yields 

= + (C6) 


where 

Jli=9ki{^'{h\)) (C7) 


Note that this population average susceptibility is the 
solution of the population averaged part of the dynamics, 
namely. 


Note that all elements of are of order 1/7V including 
the diagonal. This is becuase Xfcfcwhich is of order 1 
adds only a 1/iV contribution to the spatially averaged 
quantity, C3. 

Statistics of the full susceptibility We now turn to 
evaluate the full susceptibility,y;^^j. The random fluctua¬ 
tions in the susceptibility, i.e, the off diagnoal elements 
of X are of order l/-\/7V as will be seen below, i.e., locally 
they are much larger thany;°j. Their average second order 
moments are defined in terms of 

1 ^ 

Gkl ita,tb,tc,td) = (C9) 

i,3 

To evaluate C9, we multiply two realizations of Eq (C2) 
(differing in the time indices) and take the spatial average 
over all neurons, leading to 


Gkl i.ta,tb,tc,td) -'^Mklita - tb)Gml {taGbGcGd) = 

S {ta -tb-tc + td) 5 {ta +tb-tc- td) Ski ■ (CIO) 


where 


Fixed point At the fixed point. 


Mkiir) = gh F T)(t>'{h\{t)) . (Cll) 


Note that Gki are of order ; the contribution of to 
them is of order 1/A^ hence negligible. Furthermore, the 
uniform components of the interactions (in Eq. C2) do 
not contribute to C9 to leading order as they are smaller 
by a factor of 1/7V relative to the random contributions. 
Defining new time variables as t = ta — tb, t' = tc — td, 
T = ta + tb and T' = tc + td, equation (CIO) can be 
written as 


1 + ^1 -i]i + n{r) 


G = 2d{T-T')6{T-T')l, 
(C12) 


where ’H(t) is the P x P operator. 


n{T) = 


dr^ 


I + I-M(t), 


(C13) 


nr) = -^l + l-M (CM) 

where M is time independent PxP matrix. Fourier trans¬ 
forming C12 yields. 


[(2iD - n^)I + (w^ -h 1)/ - M] G(D, w) = 47rl, (C15) 

where D and lj are the frequencies associated with T and 
T respectively. 


G{Xl,ijj) = J dT J dTG{T,0,T,0) exTp{—iQT — iujT) 

(C16) 

. Thus, the zero frequency limit, yields a static matrix 
given (up to a constant) by 

G(D=0,w = 0) = (I-M)-\ 


and 

(see (29)). 


where M is defined in (30) and fixed point stability re¬ 
quires that all eigenvalues of M have real part less than 
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one. Note that = N = 0))^^, i.e., the square 

average of the local static susceptibility. 

Chaotic regime The largest Lyapunov exponent 
is the exponential divergent rate of the squared 
susceptibilityiV = Gfei(2t, 0,0, 0), in the 

chaotic state, which corresponds to T = 2t, and 

T' = T = t' = 0. Thus, the largest Lyapunov exponent 
is found by 


Xl = lim 4lnVGfcj(2L 0,0,0). (C17) 

t—)-oo Zt 

kl 

To evaluate it, we first write the time dependent solution 
of (C12), as 

CO 

GuiT,r, T,t') = ^r^(r, t') \r^{T)) (t')| , 

(C18) 

where {4 ’^{t)\ and |■^/)^(T)) are the left and right ^ — th 
eigenvectors of the Hamiltonian 'H(t). We note that in 
general 'H(t) is non-hermitic and the eigenvectors \tfj^{T)) 
are not necessarily orthogonal. From (C12), we find 
that the general solution for r^(T, T') is proportional to 


A^ = -i±v^r^ 

and e^is the corresponding eigenvalue of 'H(t). For com¬ 
pleteness, one must consider the boundary conditions on 
the solutions, namely that solution exists only for T > T' 
and that the second derivative of G is a delta function 
while the first derivative is finite. One obtains 


Gfci(r,T',r,T') = 

^Q[T-r)\ru{T))W{T')\ 

h 

X sinh((T-T')x/r^^ (C19) 

Finally, 


=Gfci(2t, 0,0,0) = 


4 E 


\rumwm 


e“^‘ sinh (2tyW^) . (C20) 


Thus, the maximal Lyapunov exponent is given by 

Al = —1 + \/l — coj (C21) 

where cq is the ground state of the Hamiltonian (29). 


Appendix D: Exponential transfer function 


We show that a single inhibitory population, with the 
architecture studied in Section HI and an exponential 


transfer function, (pix) = e^does not exhibit a chaotic 
phase. The fixed point of the dynamic is stable, when 
g‘^C'{x) < 1, where G = (pi), leading to 


g^ (exp (v^z + w) ) = 5^62“+^“ = g^m'^ < 1. (Dl) 

The variance, given by Aq = g^C(x), is 

^0 = 9 ^ (exp ^2 ^Aqz + 2 v}j ^ = g^m^. 

(D2) 

It follows from (Dl) and (D2) that the fixed point loses 
its stability, at the critical point gc where Aq = g^m? = 
1. For variance values higher than 1, the dynamics is 
characterized by the autocorrelation function given by 
the solution the dynamical MF differential equation in 
(45). A stable chaotic solution requires 




(D3) 


or, equivalently, the vanishing of the potential (46) at 
T = oo. Carrying the gaussian integrals in (46) one finds, 
that if the potential vanish then 

Ao = A(oo)e3^(~). (D4) 


Since every solution for the autocorrelation function re¬ 
quire 0 < A(c») < Aq, Eq (D4), can not be obeyed for a 
non vanishing variance, and no chaotic phase exist. Con¬ 
sequentially, when Aq > 1, the mean activity diverges. 


Appendix E: Dynamic equation near criticality 


Single population Using the homogeneity of the trans¬ 
fer function, the correlation function G(r) (26) can be 
written in terms of the parameter q^r) = 1 — A(t)/Ao 
giving 

C(t) = ^/l^^^(^Z + x^'^ \ . 

\ y I Z 

(El) 

Near and above criticality G(t) can be expanded in pow¬ 
ers of qir) about the fixed point value qir) = 0. Since 
linear analysis is not sufficient to account for the critical 
phenomena we keep all terms up to the first non linear 
term in (/(r), 

Aq ^“"G = {(pliz + x)) - {(pl-iiz + x))q + cqP, 

(E2) 

where c is a numerical factor of order 1. The gaussian in¬ 
tegral in diverges for 2 (z/ — n) < —1. As a result, 

the first nonlinear term may be larger than quadratic. 
The first non linear term is p = 2 for > | and p = § 
for i < < |. 
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For convenience, we redefine the gain parameter as g = 
the value of the gain g can be recovered using 
the mean field equations, (36). The dynamic equation 
(45) can be written as 


^9 = 5"(</>^(2 + a^))- 1+ 

{l-g'^v‘^{(j)l_^{z + x)))q + g'^cqP. (E3) 

Next, we expand each term in (E3) in powers of e = 
g^!g^ — 1, which controls the distance from criticality. 
We also denote 5x = x — Xc which is a function of e. Eq 
(E3) can be written as 

= a{e)+ b{€)q{T) + cqP{T), (E4) 

where 

a(e) = e + 2iygl (5a;(l + e)+ 

2v {2v — l)g^ l^z + xdl^ 5x^ + O {tSx^) (E5) 

and 

b{e) = e + 2fv%^^Hx. (E6) 

The factor is of order 1 and equals 

{v - 1) {[z + Xcf_^~^'^ v>l 

' vfc u = l . 

\{[z + Xc\+~^) + ([2 + a;c]+”^)) I < < 1 

(E7) 

In Eq (E5) through (E7) we have used the stability con¬ 
dition (22) 

v'^9l{(t>l-i{z + Xc)) = (E8) 

and the fixed point solution for the variance, which also 
holds at the transition point, (12) 

fc{(t>l{z + Xc)) = 1. (E9) 


We expect that a(e) would scale as b(e)q, implying that 
the leading order in (E5) vanish and 

Sx = -^g-^(^[z + Xc]+~^j e. (ElO) 

Thus, we can write 5(e) = eb. Following the same argu¬ 
ment as in (46), we define a classical potential for ^(t) 
by integrating the LHS of (E4), giving 


V[q] = a (e) q + l-ebq^ + 

2 P + 1 


p-ti 


(Ell) 


By applying the boundary requirements for a stable 


dV 


chaotic solution, 14[0] = V [ 9 ( 00 )] and 
gets the scaling of q and a for small e 


q{oo) = 


(1 + P) 

2p c 


= 0 one 


(E12) 


and 

a =-- —^-c( 7 '’(oo) ~ e^. (E13) 

P+1 


Multiple populations In a network composed of P pop¬ 
ulations, stability is determined by the eigenvalues of the 
PxP stability matrix, (29). Here we define a normalized 
stability matrix Mu = A~^{0)Mki = gh {cl)''^{z + Xk) 
and 


-2 _ 2 
9ki — 9ki 


MO) 

Afc(O)- 


(E14) 


We write Mki = J^fj. A^RkL* where A^are the eigenval¬ 
ues ordered according to decreasing value of their real 
part, and i?,L are the right and left eigenvectors of M 
respectively. At the transition, Ai = 1 and ReK^^i < 1. 
When gu ~ g%i, the unstable eigenvalue is Ai = 1 -f e 
where |e| <C 1 and the sign of e determines which di¬ 
rection in the space of g (away from g^j) is the unsta¬ 
ble one. To leading order, the eigenvector i?idoes not 
change. Note that Ai depends on both g and x. For 
convenience, we define e via Ai{g,Xc) i.e., as the change 
induced in Aidue to change in g for x = Xc- We are inter¬ 
ested in the properties of the system for 0 < e <C 1. For 
multiple sub-population, equation (E4) takes the form 


d'^qk 

dr'^ 


afc(e) + -f '^ghciqP, (E15) 

kl 


where = 1 — 9kfi + Xk), which vanishes at 
criticality due to (12), and Mki = gh {4’''^{z + Xk)- Since 
the matrix I — M vanishes only in the 1— direction, the 
dominant component of the vector q is in this direction. 
To see this, we make the ansatz qkir) = 1 — Afc(r)/A° = 
and assume that We then ob¬ 

tain. 


ar 2 

9r2 


ai(e) + (l-Ai)C'(r) + cV(r)^ (E16) 
a^(e) -f (1 - A^)C^(t) -f c'"C^(t)^ p ^ ^pi7) 


where a and c are the factors in the basis of M. Here 
1 —Aivanishes at the transition and is of the order of e, i5x, 
while (1 — A^) remains of order 1 near the transition. 
Thus, Q^(j) scales as C^('’')^ C{x), and the “normal 

form” of the dynamics of qk, Eq (E15), is the same as 
that of a single population, Eq (E4). In addition, in the 
generic case, we expect all qk to have nonzero projection 
on Ri. Hence, they will all scale as giving 


Qkir) 


Mr)- 


P-1 




(E18) 
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Threshold - linear transfer function In Section V, we 
have shown that in one population with threshold lin¬ 
ear transfer function ( i/ = 1), the excess mean input is 
zero at the transition (Fig. 2). Our numerical solutions 
show that in two population networks with this transfer 
function, xi are no longer exactly zero at the transition; 
nevertheless they remain small in a wide region above the 
transition. Assuming |a;/| <C 1, the stability matrix can 


be written as 

Mki = ~gliH{-xi)^Q.b~gli. (E19) 

The mean field expression for the variance, (12) can be 
written as 1 = X); gh {z + Xi))^ « 0.5 Y.i qIv Hence, 
we find that 'Yhi^ki = 1 + 0{xi), Vfc, implying that 
the eigenvector corresponding to the largest eigenvalue 
of M, is uniform to leading order m. Xk- It follows that 
near criticality, when \xk\ is small, i?^is approximately 
uniform, and all are nearly equal, as can be seen in 
Figure 10. 
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